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1 Introduction 



Effective field theories (EFTs) have become increasingly popular in describing pro- 
cesses where several scales are involved. In particular, two EFTs, namely Heavy Quark 
Effective Theory (HQET) and Non-Relativistic QCD (NRQCD) have been used for sys- 
tems with heavy quarks. These EFTs take advantage of the fact that the masses of the 
heavy quarks (charm and bottom) are much larger than the remaining dynamical scales 
in the problem. 

HQET was designed to study systems with one heavy quark |2], [| and has become 
a standard tool during the last years. Apart from the mass of the heavy quark (m) 
the remaining dynamical scales in heavy-light systems reduce to a single one Aqcd- The 
HQET Lagrangian can be organized in a power series of the inverse pole mass of the heavy 
quark. Each term in this series consists of a gauge invariant operator. Only two kinds 
of terms turn out to be important for heavy-light systems: (i) terms containing light 
degrees of freedom (gluons and light quarks) only (which are irrelevant in most of the 
phenomenological applications), and (ii) terms containing a bilinear in the heavy quark 
fields. The size of each term is easily estimated by assigning the scale Aqcd to whatever 
is not a heavy mass in the Lagrangian. 

NRQCD was designed to study systems with a heavy quark and a heavy antiquark 
[5|] and, although it is older than HQET, it has not received much attention until recently 
Pj-fOj]. In this case, apart from the heavy quark mass, there are at least two dynamical 



scales. Namely the typical relative momentum in the bound state p and the typical 
binding energy E. Because of the existence of these two scales, the power counting rules 
are different from the HQET case and the size of each term in the NRQCD Lagrangian 
is not unique. Nevertheless, counting rules have been given to estimate the leading size 
of each term (see H). Independently of the relative size of each term, the NRQCD 
Lagrangian also consists of a power series of the inverse pole mass of the heavy quark. 
Here, though, there are important terms of three kinds: the two first kinds correspond 
exactly to (i) and (ii) in HQET, where we include in (ii) terms containing a bilinear in 
the antiquark fields as well. The third kind (iii) corresponds to operators bilinear in both 
heavy quark and heavy antiquark fields (four fermion terms). 

A crucial step in building an EFT for heavy quarks is the so called matching. In 
the process of matching we enforce the effective theory to reproduce suitable S-matrix 
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elements of the full theory. In this way we fix the parameters (Wilson coefficients) of 
the effective theory. Through the matching process the high energy contributions are 
encoded in Wilson coefficients multiplying the operators in the Lagrangian (and in the 
currents) of the effective theory. The determination of some of these Wilson coefficients 
of the NRQCD Lagrangian is the main topic of this paper. 

The question arises whether the Wilson coefficients of the terms (i) and (ii) in HQET 
and NRQCD are the same. We shall support below the claim in || that this is indeed the 
case. Therefore, since the mass of the heavy quarks is (by definition) much larger than 
Aqcd, the matching may be done order by order in 1/m and a s . 

The matching for NRQCD has been known at tree level since long. This can be 
obtained by enforcing the tree level of S-matrix elements to be equal to those of QCD 
(QED) as mentioned above. For terms bilinear in the quark (antiquark) fields, this is 
equivalent to performing a Foldy-Wouthuysen transformation in the QCD Lagrangian. 
Although for HQET the matching at tree level for the bilinear terms can be carried out 
exactly as above jl2| , in most of the works it has been done somewhat differently: either 
by imposing the off-shell Green functions be equal to those of QCD M (see also [O) or 



by integrating out the 'antiparticle' degrees of freedom ||14|| . The Lagrangian obtained in 
this way is in fact different from the NRQCD Lagrangian. However both Lagrangians are 
related by local field redefinitions or by using the equations of motion || . Results for the 
matching at one loop have also been known in the HQET for some time . Nevertheless, 
attempts to perform the matching beyond tree level in NRQCD have not begun until 
recently. The main obstacle was that in NRQCD, unlike in HQET, the kinetic term was 
thought to be a necessary ingredient in the quark propagator for a matching calculation, 

1 - (i.i) 



k + ie k 2 

^0 ~ ^ H ie 

2m 

If a hard cut-off is used (fi << m), it can easily be seen that the matching can be 
performed just like in HQET since k° » k 2 /m in the ultraviolet. However, if dimensional 
regularization is used, the high energy modes {k > m) are not explicitly suppressed and 
they give non-vanishing contributions. This can be seen because the behavior of the 
NRQCD propagator changes at energies larger than the mass. In spite of this, one would 
like to use dimensional regularization because it keeps all the symmetries of QCD and, 
moreover, the calculations are technically simpler. Several authors have addressed this 
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problem |7j and recently an appealing solution has been proposed || . There, it is claimed 
that the matching in NRQCD using dimensional regularization should be performed just 
like in HQET, namely the kinetic term must be treated as a perturbation. Let us make 
some remarks which support this approach. The key point is that in order to carry out the 
matching it is not so important to know the power counting of each term in the effective 
theory as to know that the dynamical scales of the effective theory are much lower than 
the mass. The power counting tells us the relative importance between different operators 
but this does not change the value of the matching coefficients. That is, we only need 



m»\p\,E,A QCD (1.2) 

no matter what the relation between |p|, E and Aqcd is- The above becomes clear if one 
thinks of the matching as a procedure to integrate out high energy degrees of freedom 
a la Wilson: the effective Lagrangian that we obtain after integrating out energies and 
momenta until a scale /i, m » fi » |p|, E, Aqcd does not depend on the relative weight 
of the lower scales. 

In addition, in ref. dimensional regularization was used to regulate both the ul- 
traviolet (UV) and the infrared (IR) divergences in the full and the effective theory fT5 |. 
The latter arise when the S-matrix elements are expanded about the residual momentum. 
In fact, it is not so important to know the way the UV divergences of the full theory are 
regulated since the comparison is done between S-matrix elements which are UV finite 
(after renormalization) . Nevertheless, it is essential to regulate in the same way the IR 
divergences in both the full and effective theory in order for them to cancel out. This will 
always happen since by construction both theories have the same IR behavior. It is also 
important, from a practical point of view, to regulate the UV divergences of the effective 
theory using dimensional regularization. In this way, the calculation in the effective theory 
becomes trivial since there is no dimensionfull parameter in any integral. In the ref. ||, 
the matching was performed at one loop until 0(l/m 2 ) for operators bilinear in the quark 
fields. It is the aim of this paper to perform the matching at one loop until 0(l/m 2 ) for 
four-quark operators and hence to complete the matching at one loop 0(1 /m 2 ). 

We are thus faced with the computation of S-matrix elements of four heavy quarks 
in QCD and HQET. The computation of these matrix elements in HQET is unusual, 
although some related calculations already exists in the literature |T1| |I7J. Indeed, for 
heavy-light systems four fermion operators are relevant only when two of the fermions are 
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light. For heavy quarkonium systems instead all four quark fields are heavy. In fact it 
is in these S-matrix elements where we can see the peculiar IR behavior of heavy-heavy 
systems, which eventually gives rise to the Coulomb pole and hence to the standard 
non-relativistic weak coupling bound states. 

This IR behavior should appear in both the full and the effective theory. If we expand 
about the residual momentum the matrix elements of the dimensionally regulated QCD, 
we may expect an IR singularity reflecting the Coulomb pole. However this singularity 
corresponds to an odd power-like IR divergence and hence it is put to zero in dimensional 
regularization. This is not a problem. Indeed, since the effective theory has the same IR 
behavior, it also has an IR divergence reflecting the Coulomb pole which is consistently 
put to zero by dimensional regularization. The important thing when doing the matching 
is to take into account all the non-analytical behavior in the heavy quark masses which can 
not be obtained in the effective theory. Proceeding in this way we are certainly taking 
into account all the non-analytical behavior in the masses coming from high momenta 
(QCD logs). The remaining non-analytical behavior (Coulomb pole) is encoded in the 
effective theory. 

Although we have been talking about QCD and NRQCD, the results for QED and 
NRQED follow trivially from our calculations. 

Let us finally mention some of the possible applications of this work. The unequal 
mass case may be important for the B c system (this system has been studied in refs. |I8| 



which is expected to be seen in the future. This case is also important in QED for the 
muonium or Hydrogen-like atoms. For the equal mass case, our results fixes the scale of 
the a s running constant for annihilation contributions to the four quark interaction. This 
is important since in QCD, at the scales of Bottomonium and Charmonium, a s strongly 
depends on the scale. Moreover, since there are many scales in the game (m, p, E) it is not 
a priori clear which one should be used in order to fix the value of a s in the perturbative 
fl9l , |20f] and non-perturbative potentials pi] . In fact, depending on where the contribu- 



tion comes from, this value may be different. Recently, the spectrum of T(ls) and J/ip 
has been obtained from perturbative QCD at 0(maf) ||20|| . Next improvement, namely 
0(ma 5 s ) receives contributions from the matching of four quark operators, and hence our 
calculation becomes relevant. It should also be taken into account in parameterizations 



of the non-perturbative heavy quark potential along the lines of reference [21 



We distribute the paper as follows. In sec. 2 we define our four quark operators and 
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their Wilson coefficients. In sec. 3 we calculate the Wilson coefficients for the unequal 
mass case. In sec. 4 we calculate the Wilson coefficients for equal mass case. In sec. 
5 we discuss a few relevant issues in our calculation. The last section is devoted to the 
conclusions. A few technical points concerning the Coulomb singularity are relegated to 
an Appendix. 



2 Setting the matching 



The piece of the NRQCD Lagrangian containing four quark operators at 0(1/ m 2 ) 
reads 

SCnrqcd = ^1X2X2 + ^W^PlX2 (T X2 

m\m2 m\m,2 

+ ^^T> lX2 T a X 2 + -^^TV^xlTV X2 , (2.1) 
mim 2 mim,2 

where tp is the Pauli spinor field that annihilates a heavy quark and x is the Pauli spinor 

field that creates a heavy anti-quark. The subindices 1,2 denotes the possibility of working 

with different particles (different masses). We will omit these indices when the particle- 

antiparticle case is treated. 

There is another possibility of writing down these terms by using Fiertz transforma- 
tions. It reads 

d c d c 
SjCnrqcd = — —^1X2X2^1 H —^WX2X2^1 

m\VTL2 171x1712 

+ ^^t T a X2X t T<Vi + ^^t TVx2X t TV ^ (22) 
171x1712 171x1712 

The relation between the two bases is 

d^_Sd^_ Nl-l JC _ N 2 -l ^ c 
ss ~ 2N C 2N C 4iV c 2 vs 4iV c 2 vv ' 
d c d c N 2 — 1 N 2 — 1 

A - ss _|_ sv — c A c _|_ c A c 

2N C + 2N C AN 2 vs+ AN 2 vv ' 

d c 1d c 
d — -d c - 1d c + vs 4- 

d c d c 

A — _A C 4- A c J V J- - vv (0 1) 

u vv — u ss -|- u sv -|- . \^-°) 
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Of course, one can always use a redundant bases with the eight operators (2.1) and 
(2.2). The Lagrangian (2.2) is more convenient, as far as the matching calculation is 
concerned, when one is dealing with annihilation processes in the equal mass case. Nev- 
ertheless (2.1) is a better option when one addresses a bound state calculation. We shall 
use (2.1) for the unequal mass case and the redundant basis for the equal mass one in 
order to ease comparison with the actual calculations. 

In the QED case we have 

SC NRQED = —^—^1X2X2 + -^—il>i<ril>ixWx2 , (2.4) 

mim 2 77117722 

SC NRQED = — —^1X2X2^1 H —^Wx^xW^i ■ (2-5) 

77li77l 2 777x7772 

Now, the relation between the two bases is 

3 ~ 2 2 ' 

d v = -y + y- (2-6) 

We shall expand the dimensionally regulated matrix elements about zero residual 
momentum. Since there are no derivative terms in (2.1) and (2.2), the zeroth order in 
the expansion will be enough. Namely we only have to calculate the matrix element 
for the four quarks at rest. This means that the amputated legs in a digram only have 
to be multiplied either by p + (projector on the particle subspace) or p_ (projector on 



the antiparticle subspace), and the kinematic factor ym/E relating relativistic and non- 
relativistic normalizations can be put to one. We shall use MS subtraction scheme for 
both UV and IR divergences and work with the Feynman gauge. The matching coefficients 
should be gauge independent, but they depend on the subtraction scheme. It is worth 
emphasizing that we do not work in the on-shell renormalization scheme for the wave 
function (of course our masses always correspond to the pole mass), but just MS. In this 
scheme (also in MS or similars) the matching can be carried out straightforwardly. If the 
on-shell scheme is used for the full or effective theory, one must identify the UV divergences 
which correspond to a wave function renormalization and subtract them accordingly (not 
just minimally). This is obviously more tedious than using MS throughout. The little 
price to be paid for this simplicity is that our fields are not properly normalized. This 
must be taken into account by including the proper Z factors when calculating on-shell 
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Figure 1: We show the relevant diagrams for the matching of the four-fermion operators 
at order 0(1 /m 2 ) and one loop for the unequal mass case. The incoming and outcoming 
particles are on-shell and exactly at rest. 

matrix elements. 



7T V 4 v 2 W-k J ' ' 1 2N, 



c 



(2.7) 

Notice finally that at the order we are working at the Wilson coefficients in (2.1) and 
(2.2) are invariant under the local field redefinitions discussed in ||. 

3 Unequal mass case 



In this case annihilation diagrams are forbidden and we are only left with the two QCD 
diagrams of fig. 1. In these diagrams the Coulomb singularity can be identified and the 
above mentioned mechanism by which it disappears uncovered. We show this in detail in 
the Appendix. In short, things go as follow. In order to perform some integrals we have to 
move to dimensions high enough in order to regulate the IR Coulomb singularity. When 
coming back to four dimensions we can trace back the IR Coulomb singularity as a pole 
in higher dimensions, which does not appear in four dimensions since dimensional regu- 
larization loses power-like divergences. The point is that we have not provided a suitable 
dimensionfull parameter (the relative momentum) and hence dimensional regularization 
has no way to reproduce the Coulomb pole. This fact was pointed out some time ago in 

ref. urn 



We obtain the following matching coefficients 
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d„ = -C f ^-C f m? + i \-rr4 ln^ + J , (3.1) 



C/ f% " C/) ln 13 > (3-2) 



2 3 ) mf — m'2 x "~* m\ ' 



2 

mf 



mi 



m-i 



4(771? 



ms 



+ 



\ m2 1 



7^ 



«2 



1 m l 



7/- 



m 2 



777.0 



' m? 1 
ln-£ + - 
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777,1 



mim 2 



1 m 2 
In— ^ + 

v 7/2 3 



m n 



7/ z 



1 m T 

In— ^ + — 

. 7/2 3 



(3.3) 




tin ; 



+ 



2C f a 2 s , m 2 

— 5 ? mim 2 In — ? 

mf — mg m 2 



777, 



4(777? — m|) 
where Ca = N c . 

In the case of QED our results reduce to 



rrir, 



a 



777S 



2 2 

mf — 



ln^f + - 



— m 



' m? 



q 1 m i 

3m!m 2 In — ^ 

777n 



r 



(3.4) 



(3.5) 



or 



2 

mf 



mt 



9 ""1""^ — o 



r mim2 In 



(3.6) 



The spin dependent piece, which is subtraction point independent, agrees with the 
result obtained by Caswell and Lepage in ref. The scalar piece is new. 

Since we will need the equal mass results in the next section let us display them here. 
For QCD they read 



d. 



-Of 



Ca 
2 

--C 



C f )a 2 \ In 



Ca 
2 



Cf)a 



m 
2 



2 
3 



(3.7) 
(3.8) 
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d vs = -2C f a 2 s Mn ^ - ^ ) + - A C A a 2 s ( In ^ + ^ ) , (3.9) 

2 / ™2 




C A ai ( m 




d vv = 2C f a 2 s + (In ^- - lj . (3.10) 

For QED we have 

(3.11) 
(3.12) 

Recently, the scalar piece for the equal mass case in QED was calculated in ref. || 
using a cut-off regularization. This result agrees with ours except for a finite piece, which 
may be due to a different renormalization scheme for the four fermion operators in the 
effective theory. This contribution is relevant for the full calculation of the positronium 
energy levels at order 0(ma 5 ). Work in this direction is under way . 



4 Equal mass case 



For equal particles annihilation processes are allowed and they should be taken into 
account (fig. 2). From fig. 2a we obtain the well known result || 

d^ a = -na s . (4.1) 

This is the lower order (tree level) contribution. Let us consider first the one loop 
contributions arising from the gluon self-energy. Each heavy quark loop (fig. 2b) gives a 
contribution 



5« s a« , m 



2 N 



d c v f = (-7ra s )T R __•+_• In _ I (4.2) 

where Tr = 1/2 for QCD and Tr = 1 for QED. The QED result had been already obtained 
in refs. 0, ||. Light quarks (nj) and gluons give a contribution (fig. 2c) 



(4.3) 
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> 



(b) 




(c) 



(d) 




(el 





(g) 



Figure 2: We show the relevant diagrams to the matching for the four-fermion operators 
at order 0(l/m 2 ) and one loop that only appear for the equal mass case. The incoming 
and outcoming particles are on-shell and exactly at rest. 

The quark self-energy diagrams (fig. 2d) do contribute to the matching at 0(l/m 2 ) in 
the MS scheme, even though they do not in the on-shell scheme. The matching coefficient 
reads 



d 



c,2d 
vv 



, Ct a s ( 3 , m 2 
4 x (-na,)-?-— -In — - 1 



(4.4) 



2 n \4 

Let us next consider the vertex corrections (figs. 2e and 2f). Fig. 2e is quite interesting. 
A singularity associated with the Coulomb pole should appear, but again it does not show 
up when doing the computation in dimensional regularization for quarks at rest. This is 
totally analogous to what happened with diagrams in fig. 1 in the previous section. No 



signal of infinity imaginary anomalous dimension appears either [ 16| . We refer the reader 
to ref. |17j for a detailed explanation on what is going on. We obtain 

3a, ' - 2 



c,2e 
vv 



2 x 



2 s 



, m" 4 



4tt V v 2 3 
For QED fig. 2e has been already computed || |9| . 
The diagrams of fig. 2f do not appear in QED. They lead to 



(4.5) 



„9f , , / 3a s \C A f, m 2 8, 16 4 \ , A . 
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Finally we consider the contributions from the diagrams of fig. 2g. These diagrams 
also exist in QED and their contributions in this theory have already been calculated in 
ref. M. We obtain from the diagrams of fig. 2g 



= alC 



f 



-C f ) (2-21n2 + Mr) 



d c v f 9 = ^ ( - ~ + 4C f ) (2 - 2 In 2 + m 



d c ' 2g 



In 



m 



+ -(2-21n2 + «7r) 
6 



and d c ^ 9 is zero. For QED we reproduce the results in @, @|. 

Summarizing all the contributions from annihilation diagrams we obtain 



d C ss 
d° 



C A 
2 



aiC f [^-C f \(2-21n2 + i7r) 



2 



--C A + 4C/) (2-21n2 + 2vr) , 



c 

vv 



-not* 



+C A 



1 - 

11 
12 



7T V 



, m 
In — + 



1 
3 
109 



/ m , 5 
-nt n — + 2 n2 - - 
° i \ v 2 3 



m 



VK 



- + - In — 



36 



+ C f [-4] 



(4.7) 

(4.8) 
(4.9) 



(4.10) 
(4.11) 

(4.12) 



(4.13) 



Recall that we have to add to the annihilation contributions above the contributions 
(3.7)-(3.10). 

For QED we have 



d c s 
dl 



a 2 (2-21n2 + ivr) 



—Tta 



a I 44 1 , m 2 

7T \ 9 3 V 1 



(4.14) 
(4.15) 



5 Discussion 

Let us first address the important question on how the matching calculation helps to fix 
the scale of a s . In the previous section we have not paid any attention to the flavor 
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dependence of a s . For simplicity let us focus on the case of a single heavy flavor. Suppose 
that in QCD we have Nf flavors. Then in NRQCD we have rif — Nf—1 relativistic flavors. 
Consequently, the MS running coupling constant in NRQCD is expected to run according 
to Nf — 1 flavors. However, this is not obvious from the matching calculation. Notice 
that the a s s in the NRQCD Lagrangian (both explicit and in the Wilson coefficients) are 
those inherited from QCD and hence one may be tempted to make them run with Nf 
flavors. In order to clarify this issue consider first the pure gluonic part of the NRQCD 
Lagrangian ||. 



(5.1) 



where 



ds 



a* 



1 Tr In m If 

3vr 1 



-T, 



60vr 
13a s 

36(hr 



R, 



(5.2) 



T. 



R ■ 



Notice that the kinetic term does not have the standard normalization anymore. This can 
be recovered by a simple redefinition of the gluon field. Since the remaining gluon fields 
in the NRQCD Lagrangian are multiplied by g, this is equivalent to make the change 



9 -»• 9 



1 - ^T^lnm 2 /^ 2 
3n 



a* 



1 — — TRlnm 2 /^ 2 
3ix 



(5.3) 



in all the gs which multiply the gluon fields. At one loop this is nothing but changing the 
running coupling constant of Nf flavors for the running coupling constant of Nf — 1 flavors 
which is a desired result. However there are additional dependences on a s in the NRQCD 
Lagrangian (which are not multiplying gluon fields) in the Wilson coefficients. Notice 
however that the difference between a s f and af f 1 is higher order in a s and hence we 
can safely substitute a s s by a s i 1 in all Wilson coefficients except in (4.13). Indeed, the 
Wilson coefficient in (4.13) is linear in a s and hence it may be sensitive on whether this a s 
corresponds to Nf or Nf — 1 flavors. Since this a s is inherited from QCD it corresponds 
to Nf flavors. However the v dependence enters in such a way that 

Nt 



1 + 



a. 



Ti 



R 



1 m 
—nt In — 

3 i v< 



1 , m 
+ - In — 

3 z/ 2 



+ C A 



11 m 
— In —z 

12 z/2 



(5.4) 
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and hence the scale of a s is naturally fixed to m. The Wilson coefficient d c vv in (4.13) 
should better be written like 



d vv = (-7ra s (m)) 

+ c A 109 




If 5 
-nj f 2 In 2 - - - in 



36 



+ C, [-4] 



(5.5) 



Therefore, we have seen that the a s in the NQRCD Lagrangian correspond to running 
coupling constants at two different scales. The a s multiplying the gluon fields must be 
understood at some scale u, u << m and run according to Nf — 1 relativistic flavors, 
whereas the a s in the Wilson coefficients must be understood at the scale m. 

Let us next comment on the case of two different heavy flavors. If one takes the Wilson 
renormalization group point of view strictly, matching QCD with Nf flavors to NRQCD 
with Nf — 2 relativistic flavors makes sense only if mi ~ rri2. If, say, mi >> rri2, one 
should better do the matching in two steps. First one should match QCD to NRQCD 
with Nf — 1 relativistic flavors (NRQCD^-i) and next NRQCD^i to NRQCD^^- 
Nevertheless, if there is no dynamical scale between mi and rri2 and we are not interested 
in any renormalization group improvement of the Wilson coefficients, carrying out the 
matching in one step or in two steps must give exactly the same result to any fixed order 
in perturbation theory. Then, we expect our results to be useful for the B c meson in QCD 
and for the muonium and hydrogen-like atoms in QED. Recall that the Wilson coefficients 
in (5.2) trivially change into 



d 1 

d2_ 

m 2 
ck 
m 2 



d 1 = l- ^T R (Inml/u 2 + \nm 2 2 /v 2 ) , 
d 2 



d 2 

2 

m| 
d 3 



2 ' 

mi 



(5.6) 



d 2 

2 ' 2 ' 

mj wij 



Now rescaling the gluon field to its usual normalization {d\ — 1) is equivalent to 



a s —>a s (l- 7^T R (in m\j v 2 + \nm 2 2 / z/ 2 ) ^ 



-i 



(5.7) 



in the coupling constants multiplying the gluon fields, and hence these a s run with Nf — 2 



flavors. Notice also that the following equality holds, a^ 1 (is = A /mim 2 ) = as f z (^ 
v /m7m^). 



„ N f- 2 t 
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For QED an analogous discussion implies that (4.15) is given in terms of the QED 
running coupling constant. This expression in terms of the low energy a (a ~ 1/137) 
reads 

d c v = {-na) l + • (5i 




If we add (3.11) and (3.12) to (4.14) and ( |5.8p respectively, we obtain the results presented 
in 
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Imaginary parts appear in the Wilson coefficients at several instances. In order to 
obtain them from our expressions beware that we have located the cut at the negative 
real axes of the m 2 complex plain. These imaginary parts have to do with inelastic cross 
sections which cannot be obtained within the non-relativistic theory alone. They are also 
related to the decay width of heavy quarkonium states into light hadrons and had been 
calculated before [[|. Our results agree with this previous calculation. 

A word of caution is required when dealing with the Pauli matrices in D dimensions. 
The Pauli matrices arising in NRQCD have in fact very different origin, as we comment 
next. For the non-annihilation diagrams the Pauli matrices originate from 

p+a^p+Qp-a^P- =: —r]?*,^] ® K a j ] = (D — 2)a k ® a k . (5.9) 

While the first equality can be understood as a definition, for the second one we have 
used the following prescriptions (with the proper limit when D — > 4) 

[a\ a j ] = 2ie ijk a k , e ijk e ijk ' = (D - 2)5 kk ' . (5.10) 

The finite part of d vv depends on these prescriptions. For the annihilation diagrams the 
Pauli matrices originate fromQ 

p+YP- ®P-1hP+ =■ -(? k ®o k . (5.11) 

When carrying out calculations in dimensionally regulated NRQCD the same prescriptions 
have to be used and it may be eventually important to keep in mind the different origin 
of the various o k <8> o k and 1 <g> 1. For NRQED there are no ambiguities at this order since 
the spin dependent terms are finite. 



1 Notice that a k <g) a k in (5.9) and (5.11) act on different spaces even though we did not write this 
distinction between them explicitly. Recall also that the Fiertz rearrangements of section 2 only hold in 
four dimensions as well. 
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Let us finally mention that the matching at one loop for the electromagnetic current 
at leading order in 1/m arises trivially from the calculation of the diagrams in figs. 2b 
and 2c. Schematically, the full current is approximated by 

J em = _> (i + + 0{\) , (5.12) 

2 m z 

where 5 r encodes the one loop correction due to hard gluons. Now, one only has to realize 
that the relevant computation (the matching procedure follows analogously to the one for 
the four-fermion operators) is the one we performed for the diagrams above but taking 
into account the different color factors. We obtain 

* = Fb(^ + cF%^' e ) = - 40 '» 1 (5 ' 13) 

which agrees with the well known result. We consider this procedure by far the simplest 
and most efficient method to obtain 5 r (one can also trivially obtain the result for QED, 
5 r = — Aoc/tt). No problem with the Coulomb pole appears through the calculation. Notice 
also that no anomalous dimension appears either. This could be traced back to the fact 
that in both QCD and HQET for one quark and one antiquark (the effective theory to 
which we are matching to from a practical point of view) have symmetries which protect 
this current. For the effective theory this symmetry is U(A) [^31 . 



6 Conclusions 



We have calculated the matching coefficients of the four quark operators of NRQCD at 
one loop and 0(l/m 2 ). We have considered both the unequal and equal mass cases. We 
have shown explicitly how some matching coefficients in the NRQCD Lagrangian conspire 
in such a way that all a s appearing in them must be considered at the scale m whereas 
the a s multiplying the gluon fields must be considered as running with the number of 
remaining relativistic flavors only. 

The binding energies of T(ls) and J/i/j have been recently obtained from perturbative 
QCD at 0(maj) [2C|. Next improvement, namely 0(ma^) in the NRQCD framework 
requires the knowledge of the matching coefficient of the four quark operators at one loop 
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calculated here. In the framework of NRQED, these are also necessary to obtain the 
positronium binding energy at 0{ma b ). 

The unequal mass case in NRQCD may have eventual applications to the B c meson. 
For NRQED it may be relevant for precision calculations (involving recoil corrections) in 
muonium and hydrogen-like atoms. In particular it would be relevant for the spectrum 
of an hydrogen atom at 0(ma 5 ) where the electron has been substituted by a r particle. 
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A Coulomb singularity 

In this appendix we show how the Coulomb singularity is reflected in our calculation. 
Consider the following integrals 

j _ f 1 1 1 (A1) 

n J {2ti) d (q 2 + ir]) n q 2 + 2mq° + irj q 2 - 2mq° + irj ' 1 ' ; 

n = 1 and n = 2 appear in the calculation of the diagrams in fig. 2e and fig. 1 respectively. 
Upon integration over q° we obtain IR singularities from the poles of the quark propagators 
and from the poles of the gluon propagators. 

The poles of the quark propagators produce the Coulomb singularity 

' i V i 



J A \ q 2 I cfm 



where A — > is an IR cut-off. At D=4 this integral has odd power like singularities 
which are ignored by dimensional regularization. However we expect these singularities 
to shown up as poles in an odd number of dimensions. 

The poles in the gluon propagators also give rise to IR singularities. These read 

/ dD ^(^] 1 1 A ^ 2n - ( A - 3 ) 

7a W (a^m 2 



I G 
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For n—1 and n = 2 we expect a pole in D = 4 and .D = 6 respectively (an extra pole at 
D = 4 cannot be ruled out a priori for n = 2 but it will not turn up). 

The explicit result for I n below fulfills the expectations above 

T 1 f' 1 ] 71 (m 2 \ T(n-e)T(2e-2n + l) 
n (4vr) 2 \m 2 ) \Att) T(2e + 2-n) 1 

where D = 4 + 2e. Notice that the e = n singularities are of UV origin. 
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